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Abstract 

We  present  her  an  0(n)  probabilistic  algorithm  for  computing  the 
volume  of  the  union  of  n  spheres  of  possibly  different  radii.  The 
methd,  which  is  an  application  of  techniques  developed  by  [Karp,  Luby, 
83],  can  be  extended,  in  a  straightforward  manner,  to  compute  the 
volume  of  the  union  of  n  objects  (where  each  of  them  has  an  easy 
description  e.g.  boxes  of  spheres)  in  k  dimensions.  Its  time 
complexity  is  then  0(nk).  We  also  examine  the  related  problem  of 
computing  the  number  of  spheres  (or  disks,  in  the  plane)  among  a  given 
set  of  spheres,  containing  a  given  point.  For  the  case  of  n  disks  of 
the  same  radius  r,  we  can  answer  such  a  query  in  time  O(log  n  and  0(n  ) 
preprocessing  space. 

For  the  more  general  problem  of  n  spheres  of  different  radii,  we 
can  answer  such  queries  in  O(log^n)  time  and  storage  0(n  log  n), 
following  a  technique  of  [Chazelle,  83].  This  leads  to  an  0(n/n) 
expected  time  union  estimation  algorithm. 

The  probabilistic  estimation  of  the  union  follows  ideas  developed 
by  R.  Karp  and  M.  Luby  (see  [Karp,  Luby,  83]).  Some  of  our  notation  is 
heavily  affected  by  their  notation. 

We  also  show  how  to  use  the  above  methods  to  tet  if  n  spheres  have 
a  (nonzero  measure)  intersection,  in  probabilistic  time  0(n). 


■  I . . :  ':  ->'■'  I  ■ 
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1.   Introduction 

The  problem  of  estimating  the  area  of  the  union  of  many  circles  in 
the  plane  was  first  posed  by  [Shamos,  78],  It  is  straightforward  to 
solve  this  problem  by  an  0(n  )  time  deterministic  algorithm.  We  show 
here  how  to  use  Monte  Carlo  techniques  (developed  by  [Karp,  Luby,  83] 
for  estimation  of  the  failure  probability  of  an  n-component  system)  to 
get  an  0(n)  probabilistic  algorithm.  This  algorithm  can  trivially  be 
extended  to  compute  the  volume  of  the  union  of  n  spheres  (or  boxes,  or 
any  collection  of  n  fixed  description  objects)  in  k  dimensions,  in  time 
0(nk).  The  accuracy  of  the  algorithm  is  controlled  by  the  algorithm 
implementer.  Its  running  time  is  optimal  in  the  sense  that  0(n)  time 
is  needed  to  compute  the  sum  of  n  real  numbers,  under  quite  general 
models  of  computation.  An  application  of  this  algorithm  is  a 
probabilistic  0(n)  time  method  to  test  if  n  given  spheres  have  a 
nonzero  measure  intersection.  This  can  also  be  extended  to  test  for 
intersections  in  k  dimensions.  No  efficient  algorithms  for  the  problem 
of  computing  volumes  of  unions  of  objects  in  more  than  two  dimensions 
were  presented  in  the  past.  The  fastest  up  to  now  algorithm  for 
testing  if  n  spheres  intersect  was  invented  by  [Schwartz,  Hopcroft, 
Sharir,  83]  and  runs  in  time  0(n  log  n).  This  algorithm  detects  also 
zero  measure  intersections. 

The  time  complexity  of  standard  Monte  Carlo  techniques  for 
estimating  volumes  of  union  of  many  objects  depends  crucially  on 
time-efficient  solutions  of  the  cardinality  of  point  enclosure  problem 
i.e.,  given  n  objects  and  a  point  X,  find  how  many  of  the  objects 
contain  X.  Let  us  call  this  number  the  cover  of  X.  Our  0(n) 
probabilistic  volume  of  union  estimation  technique  uses  a  method  of 
estimating  covers  which  can  be  very  slow  (up  to  0(n))  but  has  the  nice 
property  that  the  slower  it  is,  the  less  is  the  number  of  random  points 
whose  cover  has  to  be  found.  Since  the  n  objects  are  given  in  advance 
and  the  covers  of  many  points  have  to  be  estimated,  one  can  use 
preprocessing  ideas  to  solve  the  cardinality  of  point  enclosure 
problems.  For  the  case  of  n  disks  of  the  same  radius,  r,  in  the  plane, 
we  can  solve  the  cover  estimation  problem  in  O(log  n)  time,  with  0(n  ) 
preprocessing  space.  This  method  however  uses  as  a  preprocessing  idea 
the   estimation  of  all  k-Voronoi  diagrams  of  the  n  centers  for  k  =  1  up 
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to  n  and  has  preprocessing  time  O(n^log  n)  which  is  prohibitive  for  use 
of  the  method  for  union  estimation.  For  the  (more  general)  case  of  n 
spheres  of  different  radii  in  3  dimensions,  we  can  answer  such  queries 
in  0(cover(X)+log  n)  time  each,  with  preprocessing  storage  0(n  log  n) 
and  also  preprocessing  time  0(n  log  n).  This  leads  to  an  0(n  log  n) 
Monte  Carlo  technique  for  estimation  of  the  area  of  the  union  of  n 
spheres,  in  case  the  maximum  cover  is  O(log  n).  We  show  how  to  get  an 
expected  0(n  /n)  Monte  Carlo  technique,  in  case  the  maximum  cover  is 
"(log^n). 


Extended  MC 
2.   Definitions  and  useful  remarks 
2.1 •  "  Notation 

Let  C,,...,C  be  n  spheres  of  radii  ri,...,r^   and   center  vector 
-»•->■  ■  ■  n 

coordinates  w ^, , . . ,v^.        Let  U  denote  the  volume  of  the  union  U  C^  . 

Let  V(C.  )  be  the  volume  of  the  sphere  C^  .  Let  Z    be   the  sum  of   the 

volumes   of   the  n  spheres.   Let  U  be  the  estimate  of  U,  produced  by  a 

randomzied  algorithm  MC  (for  Monte  Carlo).  Following  [Karp,  Luby,  83] 
we  have: 

Definition:  MC  is  called  an  (e ,5  )  algorithm  of  time  f(n),  if 
¥  e  ,6  G(0,1),  one  can  select  a  running  time  f(e,6,n)  of  MC,  such  that 
the  estimate  U  produced  by  MC  will  satisfy 

Prob  {    ^       >  e}  _<  5  (EQl) 

Definition:  If,  for  fixed  e,  5,  f(n)  =  0(n),  then  MC  is  called  a  linear 
(e ,6 )  algorithm.  In  general,  we  call  it  a  g(n)  (e ,5  )  algorithm,  if 
f(n)  =  0(g(n))  for  fixed  e  and  5 . 

2.2  Remarks  for  2  dimensions 
The  part  of  the  plane  covered  by  the  union  of  the  n  disks 
Cp..,,C^  is  partitioned  into  nonoverlapping  segments  whose  boundaries 
are  circular  arcs.  It  is  easy  to  see  that  the  number,  m,  of  these 
segments  is  0(n  ).  Since  the  area  U  is  the  sum  of  the  areas  of  these 
segments,  one  gets  a  straightforward  0(n  )  time  and  space  algorithm  for 
computing  U,  by  keeping  and  modifying  a  list  of  the  segments  created  by 
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the  first  i  disks,  for  i  =  1  to  n.  (Each  segment  is  kept  as  a  clockwise 
sequence  of  circular  arcs  surrounding  the  segment.)  An  0(n  )  algorithm 
can  be  obtained  by  keeping  just  a  list  of  circular  arcs  surrounding  the 
current  union,  of  the  first  i  disks,  for  i  =  1  to  n.  The  i+1^"  disk 
will  intersect  the  list  in  0(i)  point  and  one  has  to  update  the  list  by 
removing  arcs  which  are  interior  to  the  new  union  and  adding  the  pieces 
of  the  circumference  of  the  i+1^"  disk  which  are  not  covered  by  the 
current  union.  This  algorithm  needs  just  0(n)  space  and  outputs  the 
boundary  of  the  union  of  the  n  disks.  (The  area  can  be  found,  from  the 
boundary,  in  time  0(n),  by  the  signed  addition  of  0(n)  integrals,  one 
for  each  arc  of  the  boundary). 


2.3  The  Standard  Monte  Carlo  approach. 

Let's   assume  we   can   find  (in  time  a(n))  a  closed  surface  S  (of 

volume  V(S))  containing  the  n  spheres.   Assume  also  that  we  have  a  fast 

way  (of  time  t (n)  per  point  selecton)  of  selecting  points  within  S  in  a 

uniform  random  manner  (i.e.   the   probability   that   a  selected  point 
■  '  ■  ■ -'  ■  ■ 

falls   into  an  elementary  volume  dv,  around  a  specific  pint  Pq  ,  is  the 
same  V  Pq  in  S  plus  the  interior  of  S,   and   is   equal   to  ,-_—.• 


^     dv). 


Also,   let  B (n)   be   the   time   it   takes  to  decide  if  a  given  point  P 
belongs  to  the  union  of  the  n  spheres. 

Standard  MC 

(1)  Select  N  >_ [  ....I    -  l)  points,  uniformly  randomly.   (Note  that, 

e  ^ 
to  do  this,  we  need  to  know  just  an  upper  bound  B  on  V(S)/U). 

(2)  For  each  point  selected  in  (1),  test  if  it  belongs  to  the  union. 
Let  M  <  N  be  the  number  of  points  found  to  belong  in  the  union. 

(3)  Output  U  =  ^  •  V(S) 

Lemma  1.   The  (above)  standard  MC  method  is  an  (e ,5  )  algorithm  of  time 

0(a(n)  +  T(n)6(n)  (B-1)  J^  ] 
Proof:  See  full  paper. 
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3.   An  0(n)  (e ,5 )  algorithm 

3. 1   How  to  improve  on  the  standard  MC 

A  reasonable  value  of  a(n)  is  0(n).    (E.g.   we   can  arbitrarily 

select   a  center,   find  the  most  distant  of  the  centers  of  the  C.'s  to 

that  center,  add  to  the  distance  found  the  maximum  radius   and  draw  a 

sphere  for   S).    If  we  don't  use  any  preprocessing  ideas,  then  B  (n)  = 

0(n),  leading  to  a  time  complexity  of  0(n  •   t (n)  •   (_J i)).    If  we 

n  U     v(S) 

choose   S   to  be  the  boundary  of  the  union  U  C^  itself,  then   .  —  =  1 

but  then  the  task  of  selecting  points  from  the  interior  of  S  uniformly 

randomly   becomes  a  hard  task,  since  a  point  belonging  to  k  >  1  spheres 

cannot  be  counted  as  a  "whole"  point.   This   leads   to   the   following 

algorithm: 

Extended  MC 

(1)  Let  S  be  the  boundary  of  the  union. 

(2)  Select  N  points  (N  >_  cn/e ^ ,  c  a  constant)  as  follows:  To  select  a 
random  point  in  S,  we  select  one  C.  at  random  and  then  we  select  one 
point  P.  ,  uniformly  randomly  within  C.  . 


-1 


(3)  Compute  M  =  E   (cover(P^)) 

j=l         -J 

(4)  Output  M. 


Lemma  2.  Let  c(n)  be  the  time  to  compute  the  cover  of  a  point.  Let 
p(n)  be  the  preprocessing  time  for  this.  Then,  the  extended  MC  is  an 
(e ,6 )  algorithm  of  time  0(n  c(n)  +  p(n)). 

(Proof  in  the  full  paper).  • 

In  Section  4  we  show  how  to  gext  an  extended  MC  algorithm  of  time 
0(n  log  n),  in  cases  where  the  maximum  cover  of  any  point  is  O(log  n), 
and  an  extended  MC  algorithm  of  expected  time  0(n  /n)  else. 
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3.2   A  linear(e ,6 )  algorithm. 

We  use  here  the  following  ideas:  (a)  We  assign  greater  probability 
of  selection  to  points  belonging  to  spheres  of  bigger  volume.  (b)  We 
estimate  covers  via  a  random  test  which  takes  a  long  time  only  if  the 
cover  is  a  small  number.  This  means  that  the  spheres  do  not  overlap  a 
lot  and  we  exploit  this  knowledge  in  deciding  how  many  points  to 
select.  The  following  technique  is  a  modification  of  the  [Karp,  Luby, 
83]  techniques. 

The  linear  MC  (LMC)  method 

BEGIN 

(1)  Let  c  be  a  constant  >  20.   Let  initially  F  =  0,  number  of  trials  = 

0,  TIME  =  0. 

(2)  Compute  the  sum  E  of  the  volumes  of  the  spheres. 

v(q 


STEP  1.   Randomly  select  sphere  C^  with  probability  p^^  =  ^^ 
STEP  2.   Randomly  select  a  point  s  S  C^ 


STEP  3  f^  =  0 


repeat  until  a  C.  is  selected,  such  that  s  G  C^: 
begin 

Randomly  select  C.  with  probability  1/n 

f  :=  f+1 

test  if  s  e  C . 

TIME  :=  TIME  +  1 


end 

f 
U  :=  1  •  E 
n 


STEP  4  number  of  trials  :=  number  of  trials  +  1 
F  :=  F  +  U 
if  TIME  <  EH   then  goto  STEP  1 

STEP  5   Output  U*  =  ,      '     \         ■  , 

number. of .trials 

END  (of  algorithm  MC). 
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Note  that  LMC  stops  in  0( — •  X(n))  steps,  where  x(n)  =  time  for 
STEP  2  +  time  to  test  if  s  e  C.  . 

For  k.  dimensions,  x(n)  =  0(k.).  Hene,  LMC  runs  in  time  0(nk). 
Note  thta  LMC  could  be  used  for  objects  other  than  spheres,  which  have 
an  0(k.)  description.  It  remains  to  prove  that  LMC  is  an  (e  ,6  ) 
algorithm. 


3.3   Analysis  of  LMC 

Definition:  Let  f^  be  a  random  variable  indicating  the  value  of   f   at 
the  end  of  the  i"^^  trial  of  LMC. 

Definition.    Let   t(T)   be  the  sum  fj^  +  f2  +...+  f-p  ,  i.e.   the  sum  of 
f's  of  the  first  T  trials. 

Definition.   Let  U  (T)  be  the  value  of  U   computed  if   the  algorithm 
makes  T  trials. 

Definition.    Let   U^   be  the  estimate  of  the  area  of  the  union  for  the 
i"=^  trial. 

Definition.   Let  T(t)  be  the  number  of  trials  done  in  TIME  =  t. 

Definition.   Let  U*(t)  be  the  output  of  LMC,  given  TIME  =  t. 


Lemma  3.   For  every  i: 

nU 


(a)  The  mean  value  of  f^^  is  --— 

2 

(b)  The  mean  value  of  ff   is  <^ a 


Proof:  See  full  paper 


Lemma  4.   Let  U^  be  the  random  variable  indicating  the  value  of   U  in 
the  i^"  round.   Then  mean(U^)  =  U  and 

mean(U?)  <_  2  Z  •  U 

Proof:  See  full  paper.  • 

Theorem  1.   For  any  e  ,(S  £  (0,1) 
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Prob{  |U*(t)  -  U|  >  eU}  _<  6 


c*  n 
where   t   =  — -  ,   c  a  constant.   (See  appendix  for  sketch  of  proof  of 


Theorem  1 . ; 


3.4  Application  to  intersection  problems 
Suppose  we  want  to  check  if  the  volume  of  intersection  of  n 
spheres  is  at  least  a»E  where  a  G  (0,1)  and  I  is  the  sura  of  their 
volumes.  By  intersection  here  we  mean  overlap  of  at  least  2  spheres. 
Clearly,  this  is  equivalent  to  E  -  U  >^  a  I  i.e.  U_<Z(1  -a).  Choose 
some  £ ,5  ^  (0,1)  and  run  LMC  to  get  an  estimate  U  of  U.  Check  if 
U*  £  E  (1-a  )(l-e  )  .  If  yes,  then  with  probability  at  least  1-6  the  area 
of  intersection  is  at  least  a*E  . 

4.   The  cardinality  of  point  enclosure. 

4.1   The  case  of  n  disks  of  the  same  radius. 

We  are  given  n  disks  Cj^,...,Cj^  of  the  same  radius  r.  We  want, 
given  a  point  X  in  the  plane,  to  find  cover(X)  efficiently.  Clearly, 
cover(X)  is  the  number  of  those  of  Cp...,Cp  which  fall  into  a  circle 
of  center  X  and  radius  r.  (This  is  like  the  disk  retrieval  query  of 
[Cole, Yap, 83]  with  a  slight  modification:  Instead  of  the  list  of  the 
centers  falling  into  the  circle,  we  want  just  their  number).  We  first 
perturb  the  centers  of  C^,...,C^  inf initessimally ,  so  that  each  of  them 
has  only  1  k^^  nearest  neighbor  (for  k  =  l,...,n).  (Similar  techniques 
based  on  inf initessimal  perturbations  were  used  in  [Schwartz ,Sharir, 83] 
to  resolve  degenerate  configurations  arising  in  other  geometric 
problems).  Then,  we  do  the  following  preprocessing:  We  build  the 
k  -nearest  neighbor  diagram  (each  region  in  this  diagram  has  the  same 
k"-  nearest  neighbor)  for  k  =  l,...,n.  It  is  easy  to  see  that  the 
k'-^-nearest  neighbor,  NN,  ,  diagram  is  a  coarsening  of  the 
superposition  of  the  (k-1 )-Voronoi  and  the  k-Voronoi  diagrams  (see 
[Lee, 82]  for  definitions  and  construction).  The  k-Voronoi  diagram  has 
size  0(nk);  if  one  examines  the  method  of  constructing  the  k-Voronoi 
diagram  from  the  (k-1 )-Voronoi  diagram,  he  will  notice  that  the  edges 
of  the  two  diagrams  only  intersect  at  their  endpoints.  So,  the  size  of 
the  k'-"-nearest   neighbor  diagram  is   0(nk).    This  leads  to  a  total 
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preprocessing  time  of  O(n-^log  n)  and  total  preprocessing  storage  of 
O(n^).  To  find  cover(X),  given  X,  we  first  search  NN  ^  (i  =  1,2,...) 
until  a  2  -nearest  neighbor  of  X  is  found  at  distance  >  r.  Then,  we  do 
a  bounded  binary  search  between  NN  ,  _,  and  NN  ■  .  Let  the  furthest 
nearest  neighbor  N  of  X,  for  which  distance (N,X)  <  r,  be  the 
j'-'^-nearest  neighbor.  Then,  cover (X)  =  j.  Clearly,  the  above  procedure 
takes  O(log  n)  query  time,  per  point  X. 

4,2  The  case  of  n  spheres  in  3D  \ 

We  use  here  a  construction  of   [Chazelle, 83] .   We  are   given  n 

spheres   (of  different  radii,  in  general),  in  3D  and  we  want  to  compute 

cover(X)  for  points  X.  We  first  select  an  origin  0  and  convert   the 

spheres   into   boxes   in  polar   coordinates   r,()),9.   These  boxes  are 

parallel  to  the  axes  r,4>,9.  We  project  the  n  boxes  on  one  of  the  axes 

(say   r)   and  organize   the   projections   into  a   segment   tree   (see 

[Bentlye,80]  and  [Chazelle, 83] . _  With  this  representation,  each  node  v 

of   the   segment   tree  spans  an  interval  l(v)  given  by  the  union  of  its 

leaves'  intervals.   Let  w  be  the  father  of  v.  Let  R(v)  be   the   set   of 

boxes  whose   projections   cover   I(v)   but   not  I(w).   We  define  now  a 

"recursive"  data  structure  T2(S)  (where  S  =  {Cp..,C^})   consisting  of 

the  above  segment  ree,  and  for  each  v  a  pointer  to  the  structure  T2(V) 

where  V  is  the  projection  of  R(v)  to  the  plane  (ij)  ,9  ). 

A  query  for  cover(X)  will  now  be  answered  by  converting  X  to  polar 

coordinates   and  searching  for  the  r-coordinate  in  the  segment  tree, 

applying  the  algorithm  recursively  to  each  structure  T2(v)  encountered. 

One  can  easily  show  the  query  time  to  be  0(cover(X)  +  log  n),  with  both 

storage  and  preprocessing  time  to  be  0(n  log  n).   If  it   is   guaranteed 

that   for  all  X,  cover(X)  <  log  n,  then  the  above  procedure  gives  us  an 

extended-MC  (£,5)  algorithm  of  time  0(n  log^n). 

If,  however,  the  spheres  may  overlap  a  lot,  then  the  best  we   can 

do   is   stop   the  recursion  as  soon  as  it  seems  that  cover (X)  >^/n,  and 

run   the  probabilistic  cover  algorithm  of   LMC.    Recall   that   the 

(expected)   time  of   that   test   is  — _---  ,  leading  to  an  0(/n)  per 

_  cover(X) 

query,  and  to  an  0(n  /n)  extended  MC  (e ,6  )  algorithm. 

Further  Research 

We  pose  as  open  the  problem  of   getting  a  cover  algorithm  of 
preprocessing  time  <  "/ n   and  query  cost  <  /n  . 
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Appendlx 

Proof  sketch  of  Theorem  1  (see  full  paper  for  a  complete  proof):  We 
employ  here  the  technique  of  [Karp,Luby , 83] .  If  T  is  the  total  number 
of  trials,  then 


1   T  ,  ,  ,   -r  i-p 

i.e. 


f ^  +  . ..  +  fj  =  t(T) 

T  •  mean(f^)  =  mean(t(T))  . 


c^n 
So,  if  t  =  -i^ — _  ,  for  a  constant  ci  ,  we  expect  to  do  a  number  of  trials 

c^n    J.     c^n  r 

1  -  — -  •  7, —  <  — -  ,  because  —  <  n. 

c^Z 
Let  k  =  — — -  ("Average"  number  of  trials). 

en 
If  t"  (t"  =  — _  ,  where  c  =  C]^(l-f6))  is  the  actual  running  time  of 

6g2 
the  algorithm,  let  R(t")  be  the  event 


Clearly 


But 


|U  (t")  -  U|  >  eU 

Prob{R(t")}  =  Prob{R(t")  and  T(t")  <  k} 

+  Prob{R(t")  and  T(t")  2  k} 

Prob{R(t")  and  T(t")  <  k}  _<  Prob{T(t")  <  k} 

_<  Prob{f^  +  ...    +  f^   >   t"} 

26e  2 
_<  — ' — -  (by  Chebyshev  inequality) 

c^e^ 


Also 


Prob{R(t")  and  T(t")  >   k} 

Ui+. ..+U 
=  Prob{T(t")  2  1^  and  3  r  2  ^^^  I  "  U|  >  e  U} 

u,+. ..+u_ 

<  Prob{  3  r  >  k :  I  -i i  -  U  I  >  e  U} 
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Here  one   can  use  Kolmogorov's  Inequality  (see  e.g. 
prove  that 


[Feller, 57])   to 


,    U,+.  ..+U  gf 

Prob{    3r  >  k:        Ji -  U|    >  e  U}    <  —  . 


1 


(See   full  paper. ) 

We  conclude  that 


Prob{R(t")}  < 


26e^   86 


cy        ^1 


6  ,2e 


Let  6 '  =  __  ( — — -  +8).   So,  for  any  e ,6 '  e  (0,1)  we  can  choose  c  = 
c^(l+0  )  and  e 


such  that  if  t"  =  .'^",  then  LMC  is  an  (e  ,5  ')  algorithm. 


6'e2 
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